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Abstract. In this work we present a formulation of the Teukolsky equation for 
generic spin perturbations on the hyperboloidal and horizon penetrating foliation 
of Kerr recently proposed by Racz and Toth. An additional, spin-dependent 
rescaling of the field variable can be used to achieve stable, long-term, and accurate 
time-domain cvohitions of generic spin perturbations. As an application (and a 
severe numerical test), we investigate the late-time decays of electromagnetic and 
gravitational perturbations at the horizon and future null infinity by means of 2+1 
evolutions. As initial data wc consider four combinations of (non-)stationary and (non- 
)compact-support initial data with a pure spin-weighted spherical harmonic profile. We 
present an extensive study of late time decays of axisymmetric perturbations. We verify 
the power-law decay rates predicted analytically, together with a certain "splitting" 
behaviour of the power-law exponent. We also present results for non- axisymmetric 
perturbations. In particular, our approach allows to study the behaviour of the late 
time decays of gravitational fields for nearly extremal and extremal black holes. For 
rapid rotation we observe a very prolonged, weakly damped, quasi-normal-mode phase. 
For extremal rotation the field at future null infinity shows an oscillatory behaviour 
decaying as the inverse power of time, while at the horizon it is amplified by several 
orders of magnitude over long timescales. This behaviour can be understood in terms 
of the superradiance cavity argument. 
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1. Introduction 

The Teukolsky equation (TE) [1, 2] describes linear perturbations of a field ^' of spin s 
around the Kerr black hole solution. Black hole perturbation theory [3] has widespread 
applications in general/mathematical relativity and, ultimately, in astrophysics. To 
give a few examples, the study of TE solutions is related to fundamental topics like 
the stability of Kerr [4, 5, 6, 7], no-hair theorems [8, 9, 10], black hole quasi-normal- 
modes (QNMs) [11], and gravitational radiation from binary systems within the post- 
Newtonian approximation [12], the test-mass [13, 14, 15, 16], or the self-force (see e.g. [17] 
for a recent work and [18] for a review) approximations. 

The TE is separable if a Fourier transform in time is performed simultaneously 
with the one in the angular azimuthal direction (axisymmetric background). Exploiting 
this fact, most of the numerical calculations have been perfomed in the frequency 
domain. Time-domain numerical solutions of the TE are mainly limited to the scalar 
case (s = 0), which is extensively investigated also in multidimensional simulations, 
see e.g. [19, 20, 21, 22, 23, 24, 25]. The relevant case of gravitational perturbations 
(s = —2) has been considered for the first time by Krivan et al. [26]. There, compact- 
support initial data have been evolved with the 2+1 mode-decomposed homogeneous 
TE using a particular first order reduction and a Lax-Wcndroff second order scheme. 
The same approach has been extended to fourth order accuracy in the radial direction 
in [27], and considered in comparison to an approximate 1+1 approach in [28]. To 
date, the numerical scheme of [26] is, to our knowledge, the only successful method 
to solve the TE in the time domain for a generic spin field, and it has been applied 
in late-time decay studies of gravitational perturbations (at finite radius) [26, 29, 27], 
and in simulations of gravitational waves from a particle plunging into a rotating black 
hole [30, 31, 32, 33]. 

A long-standing problem in perturbation theory of black holes is the study of 
the late-time decay of the perturbative field, the so-called "tail". Pioneering work by 
Price [34] pointed out that the scattering of a scalar field off a non-rotating black hole 
is characterized by a power-law late-time decay oc t^^, the decay rate depending on 
the particular multipoles of the radiation considered and on the initial data employed. 
The behaviour can be related to the asymptotic form (large radii) of the black hole 
potential. Similar considerations hold for the gravitational case described by the 
Regge-Wheeler-Zerilli equations [35, 36] which have been studied in several works, 
e.g. [37, 38, 39, 40, 41, 42]. Working in null-cone coordinates, Gundlach et al. [43] 
have computed for the first time the decay rates at future null infinity (J^^) of scalar 
and electromagnetic perturbations of non-rotating black holes. It was found that they 
are different from those at finite radii. 

The numerical work of Krivan et al. [19, 26] considered for the first time the problem 
of late-time decay of fields on Kerr. As expected, the rotation of the background causes 
mode coupling (or mode mixing) between the polar modes, labelled by the integer I. In 
this paper we will often refer to the projections of a field against these polar modes as the 
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projected modes. Given an initial angular profile, specified by a particular Z'-multipole 
{pure multipole initial data), all the other projected modes allowed by symmetry/parity 
arguments can be excited during the evolution. Each projected mode decays with a 
specific rate, fii, while the main (unprojccted) field is dominated by the slowest one. 
Both s = and s = —2 perturbations were studied in [19, 26], but numerically the 
decay rates of the projected modes were studied only for the s = case. 

Tails of generic spin perturbations around rotating black holes have been studied 
analytically in two series of works by Hod [44, 45, 46] and Barack and Ori [47, 48, 10, 
49, 50]. Power law tails were studied at the horizon, finite radius and null-infinity for 
compact support, non-stationary initial data with angular dependence given by both 
pure /'-mode and generic profiles. The decay rates at the horizon, finite radius and 
null-infinity generically differ. Also, they depend on the spin-weight of the field. Mode 
coupling can be understood in the expansion of the spin-weighted spheroidal harmonics 
in terms of spin-weighted spherical harmonics. Such an expansion emerges naturally 
in a time domain analysis, whereas in the frequency domain the TE is separable and 
the spin-weighted spheroidal harmonics (depending on the frequency parameter) are 
eigenfunctions of the angular operator. Another observation related to mode coupling 
is that the decay rates of the field depend on the index I' of the initial multipole; this 
fact was referred to as the memory effect in [51]. 

Burko and Khanna [51] have reported numerical experiments in which the decay 
rates differ from the above predictions. In particular their data, computed in ingoing 
coordinates, supported a simple picture in which each projected mode decays with the 
Price law. The late-time dynamics is then dominated by the smallest of the projected 
modes which is excited (i.e. allowed by symmetry and initial data.) For s = and 
axisymmetric (m = 0) perturbations the two predictions start to differ for I' > 4. 

Further theoretical and numerical studies in the s = case [52, 20, 53, 54, 55, 56] 
helped to clarify that the decay rates depend upon the slicing of the background and 
the initial data. Hence there is no contradiction between the different predictions. In a 
recent work Burko and Khanna [56] proposed the general formula 



for decay rates of compact support, non-stationary, pure /'-multipole initial data at i+. 

Decay rates at J^~^ for s = —2 perturbations of non-rotating black holes were 
computed by Zenginoglu using hyperboloidal foliations [42] . The same technique allowed 
the first calculation of s = perturbations of rotating black holes (Kerr) at null 
infinity [57] . Recently, Racz and Toth [24] (hereafter RT) performed a detailed numerical 
analysis of s = decay rates on Kerr confirming the formula in Eq. (1) for the decay 
rates of the projected modes on the horizon and at finite radius (at J^~^ different formulas 
hold, see [44, 57, 24]). Their work is based on a particular choice of coordinates 
(see also [58, 59, 60]) that realizes a horizon penetrating, hyperboloidal foliation of 
Kerr. To our knowledge this is the only foliation of Kerr present in the literature that 




l' + l + 3 , l>l' 

v + 1 + 1 , i<r 



(1) 
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smoothly connects the horizon with future null- infinity. References [24, 25] pointed out 
an intermediate "splitting" behaviour of the decay rates of scalar fields on Kerr at finite 
radius. The phenomenon was clarified in [61] using a horizon penetrating, hyperboloidal 
foliation constructed with transmitting layers (see [62] and references therein). 

Finally, Andersson and Glampedakis [63, 29] have shown that the late time decay of 
non-axisymmetric perturbations of rapidly spinning black holes is not characterized by a 
Price-like power law but dominated by the QNM contribution. In particular, for nearly 
extremal black holes, the collective effect from several, slowly damped QNMs results in 
an oscillatory and weakly exponentially damped signal. In contrast, for extremal black 
holes the weakly damped QNMs combine to give asymptotically an overall decay rate 
oc t-^. 

In this paper we consider a novel formulation of the TE on the hyperboloidal and 
horizon-penetrating foliation of RT, and perform numerical experiments in the time 
domain studying the late time decay of generic spin perturbations of Kerr. The use of 
hyperboloidal and horizon penetrating foliations with compactification of null infinity 
has two well-known properties, e.g. [64, 65]: (i) J^"^ and the horizon are included in 
the computational domain, thus allowing an unambiguous extraction of the radiation, 
(ii) artificial boundary conditions are not needed for the solution of the Cauchy problem. 
The late time decay problem is a severe test for the robustness and accuracy of a 
numerical scheme for the TE. Note also that, for instance, Pricc-likc power law tails 
can not arise if Sommerfeld-type artificial boundary conditions arc employed [66]. Here, 
we compute numerically, for the first time to our knowledge, late-time decay rates for 
gravitational and electromagnetic perturbations of Kerr both at the horizon and at 
future null infinity. 

The paper is organized as follows. In Sec. 2 the novel formulation of the TE on the 
RT foliation is described. The key addition to the construction of RT is a s-dependent 
rescaling of the master variable, resulting in equation coefficients regular over the whole 
domain of integration. A two-parameter family of (RT-likc) coordinate transformations 
is also given. In Sec. 3 the numerical method employed is detailed. In Sec. 4 we give 
an overview of the numerical experiments performed, including convergence tests and 
a discussion of numerical difficulties. In Sec. 5 results for power law decay rates are 
presented in the axisymmetric case (m = 0). In Sec. 6 results for late time decays 
for a few non-axisymmetric cases are presented (s = 0,-2). We conclude in Sec. 7. 
In Appendix A the TE coefficients are stated in RT coordinates. In Appendix B the 
analytic calculation of the Green's function of [46] is reviewed. 

Geometric units (c — G — l) are employed. 

2. Teukolsky equation on the RT foliation 

In this section we review the construction of the RT coordinates, present a two-parameter 
generafization and derive a form of the TE with coefficients regular at and the 
horizon for every spin field value s. In the scalar case, s = 0, the equation reduces to 
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Figure 1. Conformal diagram for ttie Scliwarzscliild spacetime {a — 0). The r-slices 
(Kerr ingoing coordinates) are the dashed blue lines, the T-slices (RT coordinates) are 
the solid green lines, the t-slices (BL coordinates) are the dotted black lines. 



the one given in RT, modulo an overall factor. 
2.1. The RT foliation 

Racz and Toth [24] have proposed a hyperboloidal foliation of Kerr that penetrates 
the horizon and reaches future null infinity. The foliation can be explicitly constructed 
from the Kerr solution in Boyer-Lindquist (BL) coordinates = {t, r, 6, 0} by two 
successive coordinate transformations. To fix the notation, we recall the Kerr metric in 
BL coordinates 

5,BL = - 1 dt^ sm^ edtdip+ -^dr^ + p dO^ 

\ P J P A 

+ f + + J sin^ 9 dcp^ , (2) 

where M and aM are the mass and angular momentum of the black hole, p = 
+ a^cos6'^, and A = — 2Mr + := (r — r+)(r — r_). The two step procedure 

to construct the RT foliation is the following. First, switch to ingoing-Kerr coordinates 
= {r, r, 9, ip} defined by 

T =t-r+ / dr ^ (3) 

/dv 

Second, introduce the coordinates x"" = {T, R, 9, ip} such that 

r = T-4Mlog(|l-iJ'[) + ii|J (5) 
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The transformation in Eq. (5)-(6) realizes a hyperboloidal and horizon penetrating 
fohation with compactified Bit R — 1 [58, 59] and the r+ horizon located at 

— 0? — a? + 2Af 2 + 1 — 2 
^~ 2 (VM2 - a2 + M) ■ 



Note that R+ depends on M and a. In Fig. 1 (compare Fig. 1 of RT) the two foliations 
of the spacetime are shown in the non-rotating a = limit for simplicity. 

Following [58, 59, 24] wc can also write Eq. (5)-(6) in the more general form 

T = T + ^Jk'^ + {R/VL{R)Y -AM log(20(i?)) (8) 
R /„s 

where Vts{R) is a general compress function [65] (eventually containing a parameter S 
to fix the coordinate location of and k is a parameter determining the coordinate 
speed of the outgoing characteristic at null infinity (in RT k = S = 1). 

2.2. Scaling the TE master variable 

The homogeneous TE in BL coordinates reads 

D'^dtt^ = -A'^AaMrd^t^ + A-^9^(A"+^9^*) + sin ^-^9^ (sin ^9^*) 

+ (sin e-'^ - a^A-^)d^^^ - 2s{A-^M{a^ - r^) + (r + ia cos ^))9t* 

+ 2s(A-^a(r -M)+i cot 6 sin 9-^)8^'^ - s{s cot - 1)* , (10) 

where = (r^ + a'^)A^^ — sin^^^. The introduction of a conformal compactification 
like that of Eqs (8)- (9) requires to rescale the field variable like 

* = r-('^+^V , (11) 

in order to avoid the singularity of the physical metric at that is induced also in 
the wave equation [64, 24]. However, once Eq. (10) is expressed in RT coordinates that 
include the horizon, some coefficients have singularities on the horizon if s 7^ 0. The 
simple rescaling 

* = r-(^^+^) (Ar-^y'i/j^ A-V-V , (12) 

or 

cures the singularities. Note that the factor 1/r^ in the term (Ar"^)"* is essential since 
A ~ for r — > 00. The rescahng in Eq. (12) has been, in some form, widely used in the 
literature about the TE since [2], and thus should be considered as the "natural" one. 
In connection with the hyperboloidal transformation, it has been recently considered in 
the frequency domain framework of [67] (see Eq. (13) there). For time-domain studies 
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it is employed here for the first time together with the RT coordinates (Eq. 13). The 
resulting TE is of the form 

Cotjj + CtOtiP + CRdRijj + Cedeip + C^d^i) + CttOttiI' + CuRdRRip 

+ Ceedeei^ + C^^dy^^ip + CR^dR^i) + CtrOtr^j + Cx^dT^^Jj ^ , (14) 

where the coefficients are given explicitly in Appendix A. As one can check from direct 
inspection, the coefficients are rational polynomials and they are regular over the domain 
R G [-R+,1]. This property holds also for the general rcscaling in Eqs. (8)-(9). Note 
that the new variable ip is asymptotically constant, for example for s = —2 it is related 
to the Weyl scalar ^4 by ^ = r \E'4, since ~ ^'4 asymptotically. Note also that the 
equation in the time domain is not separable in the standard way due to the term oc d^f. 

For numerical applications Eq. (14) can be decomposed in 2+1 form separating 
each Fourier m-mode in the azimuthal direction and exploiting the axisymmetry of the 
background. The resulting 2+1 equation has the form 

Coi^ + CTdTi> + Cedei/j + CrOri/j + CrTdrTiJ + CrrOrriIj + Ceedeeil^ + CtrOtriIj = , 

(15) 

with coefficients presented in Appendix A, and the index m in the field variable ijjm 
been suppressed for brevity. 

2.3. The initial-boundary value problem 

From a PDE point of view, Eq. (14) (as Eq. (10)) is a 3+1 linear wave equation with 
variable coefficients in spherical coordinates. The principal part of the equation is 
independent of the spin weight, hence symmetric hypcrbolicity and well-posedness of 
the Cauchy problem follow from the scalar s = case when appropriate boundaries are 
considered J. 

An advantage of using a hyperboloidal foliation and compactification is that no 
artificial timelike outer boundary condition is needed, e.g. [65]. Furthermore, because 
the foliation is horizon penetrating, also the inner boundary condition is not needed. 
Inspection of the equation at i? = 1 and R = R^ gives that both boundaries are 
"outflow" boundaries. Taking the limits R — ?> R^ and i? — )■ 1 of Eq. (14) with M = 1 
(and omitting the angular terms and the non-principal part for simplicity) one obtains 

dr [pTTdr'^ + CtrOr'^^ = . (16) 

The speeds, Ci^(i?, 9) = Ctr/Ctt, for a = are given by 

CR,4R+, |) ~ -0.0460655 , (17) 

CR,+il, |) = 0.08 . 

I Numerical instabilities can still arise from unstable discretization, stiff non-principal terms, or 
exponentially growing continuum modes triggered by numerical noise. 
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Figure 2. Radial coordinate speeds for M = 1 and a = 0.9. Different panels show 
the functional dependences on R and 9, c±{R,9). Because the incoming (outgoing) 
coordinate speed vanishes at the outer (inner) boundary artificial boundary conditions 
are not needed. Note the weak dependence on 9 ioi a = 0.9. 



Similarly, one obtains for a = 0.9 

c,j,_(i?+,|) ^ -0.0723329 , (18) 

c,j,+ (l,|) ^ 0.0826788. 

The speeds cii^± are the coordinate speeds of light in radial direction that can be in 
general calculated as 

c^,+ = + Vy^« , cn,- = - V^a , (19) 

where a is the lapse, (3^ the radial component of the shift vector, and the radial 
component of the spatial 3- metric in RT coordinates (see Appendix A.) The speeds are 
illustrated in Fig. 2 for M = 1 and a = 0.9, and agree with the above estimates from 
the limiting equation. 
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3. Numerical method 

In this section we describe two different numerical strategies implemented for the 
solution of the 2+1 TE. 

The numerical algorithms rely in both cases on a method-of-lines approach for the 
time integration. Equation (15) (or Eq. (14)) is written in a first-order-in-time form, 

and evolved as an ODE system. Above -u is a certain vector of real functions and 
R{.) a discrete representation of the right-hand-side of the system. A standard fourth 
order Runge-Kutta integrator is employed. The time step is choosen according to a 
Courant-Priedrich-Lewy (CFL) condition of type At — Ccfl iiiin(/iij, /i^t), where is 
the minimum grid spacing in direction x and the factor Ccfl accounts for the maximum 
speed of the system. The spatial discretization is perfomed in two different ways, and 
the specific form of the first-order-in-time systems also differs. 

The PS-Teukode implements a fully first-order reduction of Eq. (15) with reduction 
variables^ u = {tb, drip, d^ip, dgip} and pscudo-spcctral (PS) representation of first 
derivatives. The radial direction is covered by a Gauss-Lobatto grid and Chebychev 
polynomials are employed for the PS derivatives. The angular direction 6 G (0, tt) 
is represented on a staggered equidistant grid with double covering (i.e. the function 
is extended to ^ e (0, 27r)) in order to employ the Fourier basis for the derivatives, 
e.g. [68] . Differently from [68] , here we extend the function by imposing a given parity. 
Specifically, by inspection of spin-weighted spherical harmonics, we assume the field ip 
has parity vr = (—1)™+* across the axis. We experimentally found that this prescription 
is important for long-term stability. Note that this assumption is compatible with the 
use of pure multipole initial data in the 3+1 equation and with generic initial data in the 
2+1 equation. The PS-Teukode additionally implements the following options: (i) finite 
differencing derivatives in the radial direction (as described below); (ii) the 3+1 Eq. (14) 
by using the Fourier basis also in the azimuthal direction; (iii) different fioating-point 
arithmetics: double, long-double and quadruple. Although quadruple precision requires 
several changes in the code, and since it is not available in standard hardware leads to 
a slow-down by a factor of ~ 50, it turned out to be essential for the investigation of 
s = —2 decay rates. 

The FD-Teukode implements a second-order reduction in space of Eq. (15) with 
reduction variables u = {ip, dr'ip}, and finite difference (FD) representation of the 
derivatives up to sixth order accuracy. The stencils in the radial direction are centered 
in the bulk of the domain and lop-sided at the boundaries. The angular grid is staggered 
(the same as in the PS code) and ghosts points are employed to implement the boundary 
conditions on the axis. The ghosts points are filled according to the parity condition 

^ The V variable is complex, the split into real and imaginary part is understood and omitted in the 
formulas for brevity. 
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Figure 3. Evolution of the perturbation field at the horizon and ^+ {6 = tt/2). The 
field is characterized by the quasi normal mode ringdown and a power law tail. The 
plot refers to a simulation of an axisymmetric gravitational perturbation (s = — 2 and 
m = 0) with IDO, T = 2 and a = 0.9. 

71 = (—1)™"'"'^. Artificial dissipation operators are also implemented and employed for 
long term stability. 

Stable, long-term evolutions are obtained with both codes. In this paper we mainly 
restrict attention to the results obtained with the PS code (in some cases used with FD 
in the radial direction as indicated below) since the late time decays can be computed 
more accurately with this approach. 



4. Overview of the numerical experiments 

In this section the initial data employed, the code's convergence properties, and the 
methodology used in simulation data analysis are discussed. The black hole mass is 
set to M = 1 from now on in this paper, and the spin parameter is identified with the 
angular momentum, a = J/M"^. 

Initial data are given specifying an angular profile, a radial profile, and the time 
derivative of the field. The angular profile is usually prescribed by a spin-weighted 
spherical harmonic multipole, i.e. ipi^O) oc '^Yi/^'iO) (pure multipole initial data). Four 
initial data configurations are considered corresponding to different choices of radial 
profiles and the time derivative of the field. They are named IDO, IDl, ID2, and IDS, 
and defined by 

IDO / ID2 
IDl / IDS 



>(0,i?) 
5t^(0,/?) 



G{R) I 1 




G{R) I 1 



(21) 
(22) 



with 



G(i?) = e-t(«-«o)^ 



(23) 
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These initial data are commonly used in the htcrature in which they arc referred to as 
stationary (IDO, ID2) or non-stationary (IDl, IDS), and compact (IDO, IDl) or non- 
compact (ID2, IDS) support initial data. Note that, strictly speaking, a Gaussian is not 
of compact support but we verified that, with the parameters employed, the field value 
is below the round-off level at the horizon and J^+. For example setting w = SOOO and 
Ro = 0.8, G{R) ~ 10-=^^ at the horizon (for a = 0.9 R+ 0.5222) and at i? = 1. 

We performed evolutions with the different initial data and with different spin fields 
for black hole parameters M — 1 and a — {0, 0.5, 0.9}. The quahtative outcome of all 
the simulations is shown in Fig. 3 for the exemplary case s — —2, IDO and a — 0.9. The 
figure displays the radiation field extracted at the horizon and at J^~^ for ^ = 7r/2. The 
solution has the well-known structure composed of an initial transient, the QNM phase, 
and the power law decay. We have checked frequencies and damping times of the QNMs 
in some cases 

and found discrepancies below or of the order of 1% with the tables of Berti et 
al. [11]. For example, in the case of a s = — 2 perturbation with I' — 2 on Schwarzschild 
(a = 0) the ringing frequency and the damping time are (a;, r) = (0.3876, 11.2361), to be 
compared with {uj,t) = (0.3737, 11.2410). For Z' = 3 we find (0.6133, 10.7958) in good 
agreement with (0.5994, 10.7875). For = 4 we find (0.8220, 10.6248) to be compared 
with (0.8091, 10.6202), and similarly for other cases. 

4.1. Convergence 

Let us discuss the convergence tests performed to confirm the correct implementation 
of the code. The results confirm the expected exponential convergence of the PS 
code, as far as the solution does not develop high gradients and/or reaches amplitudes 
comparable to machine accuracy. At late times the radial profile of the solution develops 
strong gradients, essentially due to the difi^erent power indices fi at the horizon and J^~^. 
The convergence rate is then slower, and progressively more coefficients are necessary 
to describe the solution. The particular shapes of these gradients depend also on s. 
Notably, in the cases s > 0, the computations are more efficiently performed with finite 
differencing in the radial direction. 

Spectral convergence is monitored by looking at the expansion coefficients, c^, of 
the spectral approximation. The highest accuracy achievable in the description of a 
smooth function is determined by round-off errors. Spectral convergence implies that the 
magnitudes of the coefficients Ck decay exponentially in k. We discuss here convergence 
of the PS code for a representative case, namely a s = —2, a = 0.9 simulation with IDO, 
I' — 2, and quadruple precision. We experimentally find that large CFL factors can 
be used without encountering stability problems, consistent with the small coordinate 
speeds of Eq.'s (17)-(18). For example here we employ Cqfl — 100, to be compared with 
the inverse of the (absolute value of the) maximum speed l/\cR^max\ ~ 1/0.0827 12.1. 
This behaviour is likely to be related to the use of the Chebychev grid (clustering at the 
extrema), the Runge-Kutta scheme, and the fact that the coordinate light speeds in the 
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Figure 4. Evolution of the spectral Fourier coefficients of 3?{^}. Given ng = 29 
angular grid points, the extended (double covering) function is represented by 2ng = 58 
complex Fourier coefficients. The real parts are the first 58, the imaginary parts the 
second 58. High-frequencies correspond to coefficients in the middle of each part. The 
left top panel shows that an initially pure multipole can be represented with a few 
points, but, as a consequence of mode mixing, higher mode coefficients are needed 
to describe the solution (top right panel). At late times only a few, low- frequency, 
coefficients are needed because only the lowest multipoles, with slower decay rates, 
have amplitudes above round-off level (bottom panels). 



bulk of the domain are very small, see Fig 2. Note that this is favorable for long-term 
evolutions and late decay studies. Figures 4 and 5 show the Fourier and Chebyshev 
expansion coefficients, respectively, of the variable The horizontal dashed line at 

-|^g-33 indicates the nominal round-off level for the quadruple precision employed. 

The plot in Fig. 4 refers to a simulation with ng = 29 angular points which amounts 
to 2n0 = 58 complex coefficients for the extended function (double covering.) The 
coefficients' real parts are the first 58, the imaginary parts the other 58. The coefficients 
corresponding to highest frequencies are placed in the middle of the real and imaginary 
sequences. The round-off level is reached in these regions, i.e. ~ 36 and k ~ 94. At 
time T = all but five coefficients are zero. This is because the spin-weighted spherical 
harmonic ~^l2o sin(6')^ = (1 — cos(2^))/2 and its derivative can be exactly described 
using the five functions {1, cos(±26'), sin(±26')}. During the evolution (e.g. T ~ 9.8, 
top right panel of Fig. 4) all other even-indexed frequency modes are excited, as a 
consequence of mode mixing. High frequency modes die off faster because they have 
longer QNM ringing phases and their subsequent power law decays are faster. That is 
why progressively fewer coefficients are necessary to describe the solution at late times 
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Figure 5. Evolution of the spectral Chebyshev coefficients of Sftj?/;} using rir = 121 
points. The top left panel illustrates that a narrow initial Gaussian needs already 
Tir = 110 points to be represented at round-off level. The top right panel shows 
that spectral convergence is achieved and at early times the solution is represented to 
round-off level. Due to the intrinsic form of the solution at late times, illustrated in 
Fig. 6, progressively more coefficients are needed to reach round-off in the tail phase 
(bottom left panel). As long as |co|/|c„^| > 10^ (bottom right panel) the solution is 
not corrupted by high frequency noise. 



(see two bottom panels). In the late stages of the simulations only the lowest modes 
survive and only very few coefficients are required. Note that due to axisymmetry 
(m = 0) only even-indexed coefficients are different from zero in 

The evolution of the Chebyshev expansion coefficients is shown in Fig. 5. The 
initial Gaussian with w = 3000 requires ~ 100 in order to be represented at 
round-off level. Every second coefficient is vanishing because the Gaussian is set to 
the center of the domain making it an entirely even function. At early times we observe 
spectral convergence, but a slower rate of convergence takes progressively over, and 
more coefficients (large k) are required to describe the function. The use of ~ 100 
does not allow to reach round-off quadruple precision at late times. This behaviour is 
due to the intrinsic form of the solution. Because the decay rate at is slower than 
at the horizon, the solution develops large gradients for i? — >■ 1, approaching a step- 
function-like form at late times. This is clearly described in Fig. 6. Experimentally we 
observe that until |co|/|cn^| > 10^ the simulations are stable, while for |co|/|c„^| < 10^ 
we observe some high-frequency noise which, eventually, corrupts the simulation and 
precludes the assessment of the decay rates. In the cases s < 0, we found that the use of 
Hr ~ 121 -=-141 points and quadruple precision allows to simulate up to T ~ 3000 -r- 4000 
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Figure 6. Late time radial profile of the solution. Due to the slower decay rate at 
J^^ the solution approaches a step-function like form at late times. This corrupts the 
spectral convergence of the Chebychev expansion. The figure refers to s = —2, a = 0.9, 
I' ~2 and IDO data. The problem is worse for the cases s > because the differences 
between the decay rates at c^"*" and finite radii are larger. 



and measure the decay rates with acceptable accuracy (see Sec. 4.2). It may be possible 
to introduce a modification of the RT coordinates that counters the "piling up" of late 
tails near the outer boundary to preserve accuracy and convergence for longer times. 

For s > resolving the radial profile of the late time solution is more challenging. 
The differences in the decay rates at J^"*" and finite radii are much larger than for s < 
cases, and a large number of spectral coefficients is required. As an example one can 
consider the decay rates of the overall field of s = ±2 perturbations with V = 2 and 
IDl. For s = —2 the rates at finite radii and J^^ are fi = 7 and /i = 6 which differ 
by only one power. The corresponding rates for s = +2 are /i = 8 at the horizon and 
= 2 at J^+, which amounts to a difference of six powers and produces large spatial 
gradients. In these situations about Ur ~ 400 points are needed in order to have stability 
at late times and measure cleanly the tail decays. Because of the quadruple precision, 
one of such simulations can take up to a month on a modern workstation (the code is 
serial). However, stable simulations at a more reasonable computational cost can be 
performed using FD in radial direction. Typical runs are performed with n,. = 801 (or 
Ur = 1601), and take about ten hours to reach T ~ 1500M, which, in turn, is sufficient 
to compute the decay rates. Obviously, also in this case, convergence is not completely 
under control and it is not possible to give a proper error estimate. However, we argue 
that these results are, at least, qualitatively correct because they are compatible with 
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Figure 7. Assessment of LPI accuracy with PS radial differentiation. The sohd blue 
line shows the difference between the LPI of the I = 1 projection at computed 
from a = 201 simulation and the analytical asymptotic value (p — —4). The dashed 
lines show the differences between the LPIs computed from simulations with different 
rir- The plot shows that the deviation from the expected value is less than 0.5% at 
late times (solid blue line). The self-differences with the target resolution Ur = 201 
are below 0.05% for = 121, of order 0.5% for = 81, and of order 5% for = 71. 
The data refer to simulations with uq = 29 of a s = — 1 perturbation with I' = 1 and 
IDl. Extraction is done at i? = 1 and 9 = tt/2. The runs were done with long-double 
precision. 



those obtained with the PS code. We compared the FD results, for several s < cases 
and few s > cases, with high resolution PS runs and found very good agreement. 
As discussed in detail in the next section, the decay rates for s > computed with 
the PS derivative and high resolution typically agree within 0.8% with those computed 
with the FD derivative. Finally, because all the s > results agree with the analytical 
expectation, we can verify their correctness a posteriori. 

4-2. Local-power-index calculation 

The power law decay of the field is monitored by the local power index (LPI) [69], 
calculated as 

^ dlog\^\ ^ ^imidT^) + ^imidr^) .24) 

^ diogT ^{tpy + ^i^y ' ^ ' 

where the reduction variable dx'ip is employed directly without taking a numerical 
derivative. Asymptotically in time, the LPI approaches (minus) the decay rate p — )• — /i. 
The LPIs of the projected modes, pi, are calculated analoguosly considering the 
projections of the fields. For a field of spin s composed of a single azimuthal m-mode. 
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Figure 8. Assessment of LPI accuracy with FD radial diiierentiation. 

The left panel compares the LPIs of the overall field of a s = +1 perturbation 
with I' ~ 1 and IDl as obtained from simulations with PS radial differentiation with 
Ur — 401 and with FD radial differentiation with = 801. Both LPIs are clearly 
approaching the analytic prediction [p = —6) (dashed blue and green lines). The 
inset shows their difference Ap = ppg — pfd ^ 10~'^ at T = lOOOM. The right panel 
compares the LPIs of the overall field of a s = — 1 perturbation with I' = 1 and IDl as 
obtained from simulations with PS radial differentiation with rir — 141 and with FD 
radial differentiation with — 801. Both LPIs are clearly approaching the analytic 
prediction (p = —5) (dashed blue and green lines). The inset shows their difference 
Ap = ppg — Pfd ^ 10^^ at T = lOOGAf. Both data refer to simulations with ng — 29. 
Extraction is done at i? = i?+ and 6 ~ 7r/2. 



ifj = 'i/'m e*'"'^, the projections are computed as 

ijiUT, R) = i'Ym^) = [ f dd sin(e) ijm{T, R, e)e''^^ %l,{e, ^) . (25) 

Jo Jo 

Note that the integration over is analytical and reduces to an overall factor 27r. The 
integration over 6 needs to be very accurate to resolve the excited higher modes / > /' 
which have small amplitudes. Using the Riemann sum or the trapezoidal rule on (0, vr) 
is too inaccurate for hq = 29 angular points. We performed a spectral integration 
in the following way. Given a periodic function on [0, 2it] and its spectral expansion 
f{6) ~ Yl,k^k€-^^^ 1 the primitive is approximated by F{9) = J d6f{9) ~ CkC^''^ / (ik) 
and the integral can be computed as F{7i) — F{0). The primitive function must be 
computed at the points 6' = 0, vr which are not included in the numerical grid, so we use 
interpolation. 

The accuracy of this technique was checked by testing the orthonormal relations 
i'^yim.l^yjm) = ^ij, which yielded the expected double precision (employed in the post- 
processing) for the projections up to / = 5 and ng = 29. 

Finally, we comment on the accuracy of the computed LPIs. Figure 7 illustrates 
how the accuracy of the LPIs is assessed. The plot refers to the / = 1 projected mode 
of a s = — 1 perturbation with /' = 1 and IDl. PS derivatives are employed in radial 
direction, and uq = 29. The solid blue line shows the difference between the LPI at J^"*" 
computed from a = 201 simulation and the asymptotic value p = —4 analytically 
predicted [46]. The deviation from the expected value is less than 0.5% at late times. 
This is a general result in all s < simulations. The dashed lines show the differences 
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between the LPIs computed from the Ur = 201 (target) simulation and others with 
different n^. The difference with = 121 is below 0.05%, with = 81 of order 0.5%, 
and with Ur = 71 of order 5%. Analogously to this case, we find that n,. ~ 121 141 
points give a reasonably accurate LPl up to T ~ 1500M for most of the simulations. 

As discussed in Sec. 4.1 for s > perturbations derivatives are computed using FD 
in the radial direction for efficiency. Specifically, sixth order stencils are employed with 
no need for artificial dissipation. The CFL factor is Cqfl = 20. To assess the accuracy 
of the LPIs, the outcomes of FD and high-resolution PS simulations are compared 
for a few test cases. A typical result is shown in Fig. 8, which refers to a s = ±1 
perturbation with /' = 1 and IDl. For s = +1 (left panel) the PS data are computed 
with Ur = 401 simulations, while the FD data with Ur = 801. In both cases the LPIs 
agree asymptotically with the predicted analytical values {p = —6) (dashed blue and 
green lines) within 0.16% at T ~ llOOM. They differ from each other by ~ 10^^ 
(~ 0.02%), as shown by the solid red line in the inset. Also for s = — 1 (and in general 
for s < 0) the FD differentiation gives accurate results, while not used because for s < 
the PS code can be used more efficiently. The right panel of Fig. 8 shows that PS data 
with Ur = 141 differ from FD data with = 801 by ~ 0.4% at T ~ llOOM (see inset). 
In this case the FD simulation looses accuracy earher and the green hne deviates from 
the blue hne. Still, the LPI is only 1.2% off the predicted asymptotic analytical value 
(p = -5) at T = llOOM. 

5. Power law tails for axisymmetric perturbations 

In this section we report the measured decay rates for s — 0,±1,±2 axisymmetric 
perturbations. In all cases we have considered pure spin weighted spherical harmonics 
multipole initial data with an axisymmetric profile, i.e. ip{d) oc *l/'m'(^) with m' — 0. 
The black hole angular momentum is fixed to a = 0.5 or a = 0.9. The parameters of 
the function G{K) are w = 3000 and Rq = 0.8. For s < simulations we used PS 
differentiation in radial direction, while for s > FD for the reason discussed in Sec. 4. 
We used grids of 141 x 29 and 801 x 29 points for the spectral and finite difference 
simulations, respectively. The CFL factor in the PS (FD) simulations is Ccfl = 100 
(C'cFL = 20 .) Quadruple (long-double) precision was employed in the s < (s > 0) 
cases. 

The decay rates of the projected modes, are given in Tabs. 1, 2 and 3 for 
s = 0, ±1, ±2 perturbations. The tables refer to the case a = 0.9. In the discussion 
of the decay rates we use the terms up/down-modes as introduced by [56]. Down- 
modes refer to the below diagonal part of the tables, i.e. / < up-modes refer to the 
upper diagonal and the diagonal part of the tables, i.e. / > The notation in these 
tables follows basically the one used by [24] but shall be clarified here. In all cases the 
asymptotic in time decay rates at the horizon and finite radii are different from those 
at future null infinity. Thus two /ii values are reported separated by a vertical line like 
x\z, where x refers to the horizon and finite radii and z to J^"*". The LPIs at finite radii 
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R+ < R < 1 vary monotonically between —x and —z at early times and reach the value 
—X only asymptotically at late times. In practice the x value in the tables is taken 
close to the horizon. In the cases s > and m = there are three distinct asymptotic 
LPls corresponding to different decay rates at the horizon, finite radii ("finite radii" in 
these cases means also "excluding the horizon"), and J^"*". The tables have thus three 
entries x\y\z. The decay rate at finite radii y is taken at a radius i?+ < R < 1 where the 
asymptotic value is reached first. Typically ^ ~ 1.05 1.2 so it is very close to the 
horizon, but the exact position is different for each data set. The larger the distance 
|it! — ^1, the greater the differences between the LPIs \p{R) —p{R)\ at early times, but, 
asymptotically in time, p{R) — )■ p{R) for any finite radius R+ < R < 1. For R > R the 
LPls vary monotonically between the values —y and —z (as long as an up-mode does 
not show splitting). For R < R the transition from y to z depends on the value of s 
(details are given in the following.) 

In some cases there are uncertainties in the assessment of the LPl due to 
inaccuracies. These can either originate from the loss of convergence in the radial 
direction due to the unfavorable form of the solution as discussed in Sec. 4.1 or to 
unresolved higher modes"'". The symbol x is employed to indicate that the loss of 
accuracy at late times prevents an unambiguous determination of the LPL In cases of 
LPI splitting we state the decay rate at finite radii in bold phase, y, to emphasize that in 
our relatively short-time simulations we can not confirm its validity for every observer. 
In particular, the asymptotic decay rate in up-modes with splitting can not always be 
assessed for observers 7? — )> 1. (In absence of splitting the LPl measurement is position 
dependent at early times, but in this case there is an unambiguous trend towards the 
same asymptotic decay rate at any finite radii.) In case of exclusion due to symmetries 
we use the horizontal line — and if we did not simulate that case we use the slash /. 

5.1. Decay rates for s — perturbations 

In this section we report the decay rates for the projected modes of s = perturbations. 
The decay rate of the unprojected field is the same as that of the lowest projected mode. 
For s = RT [24] have presented a thorough analysis of decay rates for all four initial 
data types IDO, IDl, ID2 and ID3 (and also other types). Thus for s = we restrict 
ourselves to measure decay rates for IDO and IDl to confirm the correctness of our code. 

As shown in Tab. 1 our decay rates are in perfect agreement with the work of [24] 
and of [25, 61, 56] for IDl. Figure 9 illustrates the LPl calculation for the / = and the 
/ = 2 projected modes of an /' = 4 simulation. From this example one sees that the decay 
rates of the / = projection are clearly identifiable as no = 5|4 (left panel). Instead, the 
I — 2 LPI lines bend down at radii close to the horizon and assume non-integer values. 
This is likely due to inaccuracies in higher mode projections so that in these cases we 

+ Because the decay rate /j-i* of a given projected mode I* depends on several other modes I > I* 
excited by the initial data I', the accuracy of also depends on how well the other modes I > I* are 
actually resolved. 
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ID1 , s=0, a=0.9, r=4 ID1 , s=0, a=0.9, 1=4 
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Figure 9. LPIs of the projected modes I = (left) and I = 2 (right) for a s = 
perturbation with IDl and I' — 4. Different hnes refer to different extraction radii. 
The left panel shows the decay rates approach the values fiQ — 5|4 at finite radii and 
respectively. Larger radii need longer time to reach their asymptotic decay rate 
fiQ — 5. The transition between the decay at i?+ and is monotonic. In the right 
panel there is a clear indication for the asymptotic values /i2 — 7|4 but the LPIs at 
close horizon radii depart from being a constant integer. Such cases are marked by 
brackets in the tables, here as (7)|4. 



IDl, s=0, a=0.9, r=2 
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Figure 10. Evolution of the projected modes / = 0, 2, 4 of a s = perturbation with 
IDl and /' = 2 at two extraction radii. The dashed lines refer to the modes ^ = 0, 2, 4 
at i? = 0.9410 and the solid lines at i? = 0.9883. Note that the abscissa is logarithmic. 
The plot illustrates the difference between far away observers (solid) and observers 
closer to the horizon (dashed). The observer at i? = 0.9410 measures a split of p2 at 
T lO'^ and of p4 at T 10^'^, while the observer at i? = 0.9883 does not measure the 
splitting during the simulation time. The corresponding LPIs are reported in Fig. 11. 



state the corresponding decay rates in the table with brackets, e.g. fi2 = (7)|4. 

We can unambiguously verify the LPI splitting only for the V = 2 case. Higher 
values of /' in the initial data seem to produce splitting, but we can not assess 
unambiguously the values. In agreement with previous results we do not observe 
splitting for the two lowest initial modes /' = 0, 1. The LPI splitting is described by 
Fig. 10 and Fig. 11, which refer to /' = 2 with IDl. The splitting is clearly identifiable 
by looking at log — log plots. In Fig. 10 the splitting for close horizon observers is visible 
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Figure 11. Splitting of the LPIs in the up-mode I = 2 tor s — with IDl and I' — 2 
(corresponding to figure fO). The different lines refer to 10 different extraction radii 
from R = to i? = 1. The decay rates at the horizon and are fj,2 = 5|4. At 
radii in between the LPl splits, e.g. at i? = 0.9410 (purple dashed line) the LPI is not 
a constant integer during the evolution time. Only asymptotically in time the LPI at 
any finite radius will approach the value p2 = —5. For the observer at i? ^ 0.9883 
(dashed orange line) it is possible to measure an intermediate decay rate p2 = —7. 
Such cases are marked in bold face in the tables, indicating that at intermediate times 
different decay rates than the stated ones can be measured by observers far away from 
the horizon. 

as cusps in the tail phases of the / = 2 and / = 4 modes (dashed blue and cyan lines). 
Observers far away from the horizon (solid lines) will eventually develop the cusp but 
at later simulation times. In Fig. 11 the corresponding LPIs are shown for different 
extraction radii. The asymptotic decay rates at the horizon and J^~^ are, respectively, 
/i2 = 5 and /i2 = 4, as indicated by the solid red and solid black lines. The LPIs at 
finite radii close to (dashed lines in the upper part of the plot) approach the value 
P2 = —5 the faster the closer they are to R+. On the contrary, LPIs at finite radii 
R > 0.98 have values p2 — —7 (see the black, grey and orange dashed lines) and are 
expected to approach the value at R+ only asymptotically. For these radii the simulated 
time is not sufficient to observe the splitting towards their asymptotic value. The cusp 
in the tail of the I = 2 mode at i? = 0.9410 in Fig. 10 corresponds to the dashed purple 
line in Fig. 11 which bends down at early times and "falls" from the top of the plot 
at later times. The splitting behaviour basically corresponds to a sign change in the 
field, which could be probably explained with some analytical argument. Note also that 
changing the location of the initial Gaussian pulse has a strong effect on the splitting 
behavior as discussed in [61] (see the latter reference for more insights). 

5.2. Decay rates for |s| = 1 perturbations 

In this section we present our results for s = ±1 perturbations. Table 2 summarizes the 
decay rates of the projected modes. 

Similarly to the s = case, most of the LPIs can be calculated very accurately. 
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Table 1. Decay rates /i; for s = with IDO and IDl at finite radii|futurc null infinity. 
Brackets point to uncertainties in the LPI assessment due to possible inaccuracies or 
not verifiable splitting, x to ambiguous or immeasurable values, — to modes excluded 
by symmetry and / to not simulated cases. Bold values denote splitting in time, i.e at 
intermediate times pi ^ — /i; for R < JJ^ . 
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ID1, s=-1, a=0.9, r=4 



RH=0.5222 - 
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R2=0.9605 



R3=0.9715 R6=0.9962 
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R5=0.9914 R8=0.9990 - 
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Figure 12. LPIs of the projected modes I = 2 (left panel) and I = 5 (right panel) of 
an initial s = — 1 perturbation with I' = 4 and IDl. Different lines refer to different 
extraction radii. The left panel shows the decay rates /i2 = 7|5 at finite radii and ^+ 
respectively. The transition is monotonic. The right panel shows that LPIs for the 
I = 5 projection split in time. At finite radii including the horizon and the decay 
rates are /zs = 10|8, respectively. Only asymptotically in time far out observers like 
at i? = 0.9962 (dashed black line) or at i? = 0.9883 (dashed grey line) will observe the 
LPI lines bending down and approaching ps = — 10 from above. Such splitting in time 
is indicated by bold face in the tables. 



Looking for example at the left panel of Fig. 12, which refers to the / = 2 projection of 
a s = — 1 perturbation with IDl and I' = 4, one can see that the decay rates H2 = 7|5 
can be unambiguously inferred from the plot. Table 2 reveals certain general patterns. 
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Figure 13. LPIs of the projected mode Z = 2 of a s = +1 perturbation witii V = 2 
and IDO measured by observers close to tiie horizon (left panel) and far away observers 
(right panel). From both panels one can see the splitting of LPIs into three distinct 
asymptotic values. The decay rates at R+, finite radii excluding and ^+ are 
^JL2 = 7|6|3. This splitting in space is specific to s > 0, m = perturbations. The 
left panel illustrates the behaviour of the LPIs when going from to R = 0.5347. 
Observers close to i?+ measure initially decay rates close to p2 = ^7 but asymptotically 
they will see the LPI bend down and approach the finite radii rate P2 — from above 
(in this case observed first at i? « 0.6010). Thus the splitting in space results also in a 
splitting in time for close horizon observers of s = +1, m = perturbations. The right 
panel illustrates the splitting in time for far away observers already shown for s = 
and s = — 1 perturbations in Fig. 11 and the right panel of Fig. 12. 



For example the /i/ at finite radii are always increasing* by one if one steps from one 
column to the next, i.e. to the higher projected mode. Thus the overall field at finite 
radii is always dominated by the lowest excitable mode, the / = 1 mode for s = |1|. 
This is not true for the decay at J^^, where higher projected modes can have the same 
decay rates as the / = 1 mode and can thus dominate the asymptotic signal. Neglecting 
for a moment the two lowest rows, /' = 1, 2, /x; also increases by one if one steps one row 
down. 

Another interesting result is about the effect of different initial data. Switching 
from compact to non-compact support ID results in a slower decay. To see this one can 
compare the tables for IDO with the tables for ID2 and the tables for IDl with the tables 
for IDS. Changing from stationary to non-stationary ID leaves the results for the lowest 
allowed initial mode Iq (the first row in the tables) unchanged. However, the second 
and the third rows seem to be interchanged. From the third row on the difference for 
finite radii decay rates seems, again, to be just a decrement by one when switching from 
non-stationary to stationary ID. 

There are no other numerical investigations of projected modes' decay rates for 
s 7^ 0, a = 0.9 including the horizon and J^~^, but we can compare our results to the 
theoretical predictions of [46]. In particular, the Green's function calculation provides 
a result immediately comparable with our IDl case (see also Appendix B). The decay 
rates of the overall field as stated in Tab. 1 of [46] can be compared with the first columns 
of our tables for IDl. Also, though not stated explicitly, the decay rates of projected 

* Note that increasing values of in the tables means actually the LPIs pi get more negative, i.e. the 
decay is faster. 
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modes for IDl can be read off the effective Green's functions in Eq. (15) (finite radii) 
and Eqs. (35-36) {J^'^) of [46]. In summary, as shown in the tables, our results for 
IDl are in agreement with the predictions by [46] in all tested configurations, except 
that an apparent difference can be observed for the second lowest allowed initial mode 
I' = lo + 1. Comparing with Hod's results, the values in this row are larger by one than 
the predicted values. However, this discrepancy can be traced back to a missing case 
distinction in the final parametrization of the result in [46], see Appendix B. With this 
small correction, there is complete consistency between the analytical and numerical 
results. 

The splitting of LPIs recently reported for s = perturbations [24, 25, 61] and 
verified in Sec. 5.1, is also present in s 7^ cases. It is convenient to refer to this as a 
splitting in time since intermediate decay rates pi 7^ —fii can be measured at early times 
while for T — )■ 00 a unique asymptotic decay rate /x; is obtained at any finite radii. An 
example of splitting in time is reported in the right panel of Fig. 12. The qualitative 
behaviour is equivalent to the s = case previously described, e.g. the far away observer 
at R — 0.9962 (dashed black line) measures at early times an intermediate LPIp5 —11 
different from the asymptotic value /i^ — —10. Only observers at i?+ (solid red line) and 
close to it (dashed red and green lines) measure the rate ps = — 10 during the simulation 
time. Such splitting in time is only present in the up-modes: for instance the left panel 
of Fig. 12 illustrates the LPIs of the down-mode / = 2. No splitting is visible and the 
LPIs vary monotonically between the vahics p2 — —5 (solid black line) and p2 = ~7 
(solid red line) at J"^ and -R+, respectively. 

The splitting in time differs from the splitting in the case of s > and m = 
predicted by [46, 49]. In this case the LPI has three distinct asymptotic values at the 
horizon, at finite radii, and at J^"*", which could be referred to as a splitting in space. 
The splitting in space is observed in any projected mode as well as in the overall field 
while the splitting in time is only observed in the up-modes. 

The right panel of Figure 13 shows both splittings of LPIs in the projected / = 2 
mode of a s = +1 perturbation with I' = 2 and IDG. The decay rates at and J^^ 
are, respectively, /i2 = 7 and fi2 = 3, as read off from the solid red and black lines. 
The splitting in space is clearly visibile at the radius R ^ 0.8083 which is approaching 
the value /i2 = 6 (dashed red line). The left panel of Figure 13 corresponds to the 
same simulation as the right panel but focuses on observers very close to the horizon. 
The latter measure at intermediate times LPIs close to the horizon decay rate P2 — 7 
(dashed red and green lines). Instead observers well-off the horizon measure already at 
early times the LPI p2 = —6. Note that the transition between both takes place in a very 
small spatial interval because already observers at i? = 0.5288 measure approximately 
the asymptotic decay rate p2 ~ — /i2 — —6 (dashed grey line). 

5. 3. Decay rates for | s | = 2 perturbations 

In this section we present our results for s — ±2 perturbations. 
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Table 2. Decay rates /x; for s = — 1 (left) at finite radii|null infinity and for s = +1 

(right) at the horizon[ finite radiijnull infinity. Brackets point to uncertainties in the 
LPI assessment due to possible inaccuracies or not verifiable splitting, x to ambiguous 
or immeasurable values. Bold values denote splitting in time, i.e. at intermediate times 
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Qualitatively, s = —2 simulations are in complete agreement with s = and s = —1 
simulations. As an example, Fig. 14 shows the LPIs of the / = 2 projection for I' = 4 and 
IDO at different extraction radii. Prom the plot we can extract the decay rates /i2 = 8|7, 
confirming that also for s = — 2 there are two distinct asymptotic decay rates at any 
finite radii including the horizon and at future null infinity. The radius dependence is 
monotonic as already described for the down-modes of \s\ = 0, 1 perturbations. 

The LPI calculations for |s| = 2 are summarized in Tab. 3. The decay rates 



Numerical solution of the 2+1 Teukolsky equation, application to late-time decays 25 



IDO s=-2, a=0.9, l'=4 

-6.5 n • • • 



Q. -f.O 




600 800 1000 1200 1400 1600 

T/M 

RH=0.5222 R2=0.9114 R4=0.9740 R6=0.9919 R8=0.9987 RS=1.0000 

R1 =0.8229 R3=0.9544 R5=0.9883 R7=0.9971 R9=0.9997 



Figure 14. LPIs of the projected mode Z = 2 of a s = —2 perturbation witli V = A 
and IDO. Different lines refer to different extraction radii. The asymptotic decay rates 
extracted from this plot are /i2 = 8|7 corresponding to the decay rates at finite radii 
including the horizon and future null infinity. The transition of the LPIs from to 
is smooth as for s = 0,-1, +2 perturbations. 



Near Horizon Observation - IDO, s=2, a=0.9, l'=3, m=0 
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Figure 15. LPIs of the projected mode Z = 3 of a s = +2 perturbation with Z' = 3 
and IDO measured by close to the horizon observers (left panel) and far away observers 
(right panel). From both panels one can see the splitting of LPIs into three distinct 
asymptotic decay rates at _R_|_, finite radii excluding and here /is — 9|8|3. 
This splitting in space is special to s > 0, rri = perturbations. The left panel 
illustrates the behaviour of the LPIs when going from _R+ to i? w 0.5503. Differently 
from s = +1, TO = perturbations (Fig. 13) no splitting in time is observed. The right 
panel illustrates the splitting in time for far away observers already shown for s — 
and \s\ = 1 perturbations in Fig. 11 and the right panels of Figs. 12, 13. 



follow exactly the same patterns as described in Sec. 5.2 for |s| = 1. At finite radii 
fii increases by one when switching to the next higher projection I or the next higher 
initial mode /' (neglecting again the first two rows I' = Iq, Iq + 1). The overall field 
is asymptotically dominated by the lowest mode, i.e. for m = the / = 2 projection. 
Instead, at J^"*" the / = 2 mode is not always the dominant one, because higher modes 
I > 2 can have the same decay rate. Changing from compact to non-compact initial data 
results in the decay rates decrease by one. Changing from stationary to non-stationary 
initial data leads to the same modifications of the decay rates as discussed for the case 
|s| = 1. Comparing with the predictions of [46] for IDl again all our results agree. See 
again Appendix B for the case /' = /q + 1. 

The splitting in time of the LPIs of up-modes is also observed for |s| = 2. As an 
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example the right panel of Fig. 15 illustrates the splitting in time for the / = 3 projection 
of a s = +2 perturbation with /' = 3 and IDO. The left panel of Fig. 15 refers to the 
same simulation as the right panel but concentrates on the LPIs at radii very close to 
the horizon. The LPIs at i?+ and J^~^ are, respectively, p2 = —9 and p2 = —3, and they 
are marked by the solid red and black lines. The splitting in space is clearly visible as 
the LPI is p2 — —8 at the radius R ^ 0.8095 (dashed red line in right panel). It is 
interesting to investigate the close horizon behaviour in the left panel. Differently from 
the s — +1, m — case (compare with Fig. 13, left panel), splitting in time is no^ present 
for LPIs measured by close horizon observers R+ < R. This is not a special feature of 
the particular case shown here. We observe this monotonic transition behaviour in all 
the projections, as well as the overall field, for all s = +2, m = perturbations. 



Table 3. Decay rates /i; for s = —2 (left) at finite radii|null infinity and for s = +2 
(right) at the horizon [finite radiijnull infinity. Brackets point to uncertainties in the 
LPI assessment due to possible inaccuracies or not verifiable splitting, x to ambiguous 
or immeasurable values. Bold values denote splitting in time, i.e. at intermediate times 
Pi -1^1 for R<1. 
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Table 4. Decay rates fii for s = and m = 1 with IDl at finite radii] null infinity. 

Brackets point to uncertainties in the LPI assessment due to possible inaccuracies or 
not verifiable splitting, x to ambiguous or immeasurable values, — to modes excluded 
by symmetry. Bold values denote splitting in time, i.e at intermediate times pi ^ —j^i 
for ii < 1. Square brackets point out diff'erent values compared to our m = tables. 
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Table 5. Decay rates m for s = — 2 and m = 2 with IDl at finite radii] null infinity. 

Bold values denote splitting in time, i.e at intermediate times pi =/= —fii for < 1. 
Square brackets point out different values compared to our m = tables. 
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6. Late time decay of non-ctxisymmetric perturbations 

In this section we report numerical experiments for the late time decay of non- 
axisymmetric perturbations. We consider s = and s = — 2 perturbations and 
a e [0.9, 1]. We use again pure multipole IDl, ijj{9) oc ^Y/'m', with m' — 2 [m! — 1) for 
s — —2 [s — Q) and V — 1, 2, 3, 4, 5. For s > and a — 0.9 we discuss the cases s — +1, 
V — 2 and s — +2, Z' = 3, but no further systematic investigations are performed. 

6.1. Late time decay for a — 0.9 

Let us first discuss the case a — 0.9. The numerical settings are the same as those for 
axisymmetric perturbations, but using Rq 0.76. The field dynamics at late times is 
characterized by power law tails, as in the axisymmetric case. 

The decay rates for s = are summarized in Tab. 4. Similarly, the decay rates 
for s = —2 are summarized in Tab. 5. The values different from the axisymmetric case 
are explicitly indicated by square brackets (compare with Tab. 1, 3.) In the s = case 
the decay of the overall m = 1 field can be different from the axisymmetric one, just 
because the lowest allowed mode for even I' perturbations with m = 1 is Z = 2 instead 
of Z = 0. In the s = — 2 case instead the lowest allowed mode is Zq = 2 for m = 2 as 
well as for m = 0. Still the overall m — 2 field is found to decay differently from the 
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T=4027.1 T=4027.1 




Figure 16. Spectral Fourier (left) and Chebyshev (right) coefficients of the field 
variable for s = -2 and a = 0.9999 at T « 4027M. The use of ng = 21 and 

Ur — 281 points allows to achieve spectral convergence at late times in simulations of 
perturbations of almost extremal Kerr black holes. Long-double precision is employed 
here. 



axisymmetric analogue for /' = 3. As indicated by the tables, the decay rates that differ 
from the axisymmetric case are only those relative to /' = 2 (/' = 3) for s = (s = —2) 
(second rows in m 7^ tables). See Appendix B for the mathematical origin of this. 

The s = 0, m = 1 results agree with e.g. [24] and they are consistent with various 
values reported in the literature [19, 20, 55, 56, 25], where comparable. The s = —2 
results agree completely with the formulas extracted from the effective Green function 
computed in [46]. Note that for m 7^ the parameterization in Eq. (15) and Eqs. (35-36) 
of [46] yields the correct results also for I' = Iq + 1 (compare Sec. 5). 

Furthermore, for s > the two simulations with s = +1, V = 2 and s = +2, /' = 3, 
agree with Hod's formula. For instance, in case of a s = 2 and m = 2 perturbation we 
measure fi2 = 7|2 as opposed to fi2 = 9|8|3 for m = 0, which confirms the prediction 
that the splitting of LPIs in space is not present for m 7^ 0. 

For what concerns splitting, we confirm that the splitting of LPIs in the / = 2, 4 
projections of a s = and /' = 2 perturbation observed for m = disappears in the 
non-axisymmetric case m = 1 [24]. By contrast, we observe splitting of LPIs in time in 
the / = 3, 5 projections of a /' = 3 perturbation, similarly to the m = case. 

6.2. Late time decay for a — )■ 1 

In this section we discuss the late time decay behaviour as the angular momentum of 
the background approaches the extremal value a = 1. We consider perturbations with 
s = and s = —2 and the Gaussian in the initial data has parameters Rq = 0.9 and 
w = 6000. 

Varying the angular momentum a in the range {0.9, 0.99, 0.999, 0.9999, 1} in 
otherwise identical simulations, we observe that the field develops an oscillatory 
behaviour in space, mostly localized near the horizon and of progressively larger 
amplitude for a — t- 1. The signal extracted at i?+ or J^~^ is characterized by a QNM 
phase which becomes longer (damping time increases); for a = 0.9999 a power law tail 
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Figure 17. Late time decay of s = —2 non-axisymmetric perturbations for higlily 
spinning backgrounds. The plots show the field variable extracted at i? = 

and i? = 1 and at 6* = 7r/2 for a £ {0.9,0.9999} and = 1.1345 for a = 1. Initial 
data are IDl and I' = m = 2. From left to right the field decays with a power law 
tail (a = 0.9), or with an oscillatory behaviour damped by either a slow exponential 
(a = 0.9999) or a power law l/T {a = 1). 

is not observed up to T ~ 8000 but may eventually arise at later times. This behaviour 
is consistent with the findings of [63, 29], and can be qualitatively understood in terms 
of modes trapped in the superradiance resonant cavity. 

For a — )■ 1 the amphtude of the oscillating field grows by a few orders of magnitude 
near at the beginning of the simulation (especially for s = —2). This leads to a 
step-like function which corrupts convergence if the number of radial grid points rir is 
too small (see analogous discussion in Sec. 3). Accurate simulations require PS radial 
derivatives and n,. ~ 281 points for non-extreme cases. For the extreme case we set 
Ur = 561 for s = and = 701 for s = —2. On the other hand the large amplitudes 
and small damping times allow for the use of long-double instead of quadruple (except 
for a = 0.9 and s = —2), which alleviates the computational costs. Due to these facts 
we observe spectral convergence during the whole simulation time, as shown in Fig. 16 
for s = —2 and a = 0.9999. Both the Fourier coefficients (left panel) and the Chebyshev 
coefficients (right panel) reach round-off level in approximately straight lines at very late 
times T ~ 4000. Note that the Fourier coefficients contain both even- and odd- indexed 
frequencies in contrast to the m = simulations shown in Fig. 4. 

The simulations' outcome is illustrated in Fig. 17 for s = —2, /' = m = 2, 
where we report the field extracted at the horizon and J^~^. On the simulated 
timescale, the power law tail (left panel) is observed for values a < 0.999. For larger 
values (a = 0.9999, central panel) an oscillatory and exponentially damped behaviour 
dominates the dynamics both at the horizon and J^+. For a = 1 the field remains 
oscillatory but does not decay exponentially during the simulated time. At the horizon 
the field amplitude is still growing at T = 2500, while at J^~^ the decay departs from 
an exponential law. A similar qualitative behaviour is observed in the case s = 0, 
V = m = 1. 

It is instructive to compare the QNM (complex) frequencies uj = {ur,uji) extracted 
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from the s = —2 simulations with the data by Berti et al. [11] (hsted there up to 
a — 0.9999). The fundamental (longest lived) QNM has the frequency extracted by 
counting the periods of the signal for long times. (Note that for a — )■ 1 the frequencies 
of the overtones accumulate around this value.) Damping times (or Ui) are extracted 
by fitting the envelope (maxima) of the field. For a = 0.9 we find cu ~ (0.6827, 0.0649) 
both at R+ and J^"*" which agrees very well with the target values (0.6715, 0.0649). For 
a = 0.99 we find u ~ (0.8747, 0.0294) compared to (0.8709, 0.0294) and for a = 0.999 
we get u ~ (0.9584,0.0105) compared to (0.9558,0.0105). For a = 0.9999 we obtain 
u! ~ (0.9862,0.0035), again in agreement with (0.9857,0.0035). These results also agree 
with the analytical formulas for QNM frequencies of near-extreme black holes proposed 
in [70, 71]; the discrepancy is (- 1%, ~ 5%) for a = 0.99 and (~ 0.02%, ~ 0.5%) for 
a — 0.9999. The envelope of the field for a = 0.9999 is well represented by the fit 

logio l^-^l = -1-079(3) - 0.0015064(1)T , (26) 

where the error on the last significant digit of the coefficients is reported in parentheses. 
Note that the fitting coefficient proportional to T needs to be divided by \ogiQ{e) to 
obtain Ui. In the extreme case we measure Ur ~ 0.9970 at and Ur ~ 1.0015 at 
It is interesting to notice that this is very close to the superadiance frequency 
UJ+ — ma/ {2 Mr j^) (similarly for s = we obtain Uj. ~ 0.5011 both at i?+ and J^"*" in 
the extreme case) . The envelope of the field is best fitted by the power law 

logio l^f^V^I = 1-16(1) - 1.05(2) logioT , (27) 
which confirms the expected behaviour for s = — 2 [29]. 

7. Conclusion 

In this work we have presented a novel approach to time domain numerical solutions 
of the Teukolsky equation (TE) for generic spin perturbations. The approach is based 
on the use of the hyperboloidal slicing of the Kerr spacetime, previously introduced 
by Racz and Toth (RT) [24], together with a proper rescaling of the field variable that 
leads to a regular equation for generic spin perturbations (Eq. (12)). The RT coordinates 
are both horizon penetrating and compactify null infinity {J^^) at a finite coordinate 
radius. A generalization is given by Eqs. (8)- (9), which, in particular, allows for a 
generic compactification. The TE obtained from this 2-parameter family of coordinate 
transformations and the field rescaling is regular. 

Accurate time-domain numerical solutions for generic spin perturbations can be 
performed employing standard discretization techniques. In particular, long-term, 
stable 2-1-1 evolutions are obtained with the method-of-line scheme (e.g. Runge-Kutta 
time integrators) and either high-order finite differencing or pseudo-spectral spatial 
derivatives. No boundary conditions are needed, and the computational domain includes 
the horizon and J''^. 
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As an application, we performed numerical experiments to investigate the late time 
decays of generic spin perturbations. This problem is a severe test for the robustness of 
the approach. 

In a first series of experiments we investigated power law decays of axisymmetric 
(m = 0) perturbations of spins s = 0, ±1, ±2. A variety of initial data have been used: 
stationary (IDO, ID2) or non-stationary (IDl, ID3), and compact (IDO, IDl) or non- 
compact (ID2, IDS) support. In all the cases we used pure multipole initial data. The 
numerical results for s = agree with those reported in the literature. In particular, for 
non-stationary and compact-support initial data (IDl) the decay rates follow Eq. (1). 
For what concerns the decay rates of the projected modes of s 7^ perturbations, it is 
straightforward to write down empirical formulas summarizing our results for m — 0. 
For instance we find for stationary and compact-support initial data (IDO) 

f l' + l + 3 + 5, I'^lo 
^'-\l' + l + 2 + S, l'>lo, ^^^^ 

at Above, /q = max(|m|, |s|) = \s\ and 5 = in general, but 5 = 1 in the special 
cases m = and s > 0. The corresponding decay rates at i^ are given by the same 
expressions with 5 = in all cases. At J^~^ the rates follow 

r i-s + 2 , i>i' 

^^ = \l'-s + l , Kl'. ^^'^ 
Similarly for non-stationary, compact-support initial data (IDl) we obtain 

_ / I' + 1 + 3 + 6 , l' = kJo + l 
^^-\l' + l + l + 5, r>/o + l, ^ ' 

at R^. The corresponding decay rates at i~^ are given by the same expressions with 
5 = in all cases, while at we get 

l-s + 3 , / + 1 > /' = /o + 1 = / + 1 
Hi= ^ l-s + 2 , l + l>l',l'^lo + l = l + l (31) 
/' - s , l + Kl' 

Decay rates for ID2 (ID3) can be easily deduced from those for IDO (IDl) by decreasing 
the III by one 

IDO ^ ID2, IDl ^ ID3 -.fxi^ni-l. (32) 

These empirical formulas describe all our numerical results for s 7^ and m — 0. 

The decay rates computed for s ^ can be compared with Hod's analytic 
predictions [44, 46]. The comparison is easiest to make for IDl. The results agree 
with the analytical computation (see Appendix B for a review of the calculation and a 
small correction for the case I' = Iq + 1 with m = 0). For similar type of initial data, 
Krivan et al [26] reported for s = —2 perturbations with the concerned I' = 3 initial 
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mode fi = 7.72. This value was interpreted to asympotically reach n = 7 but could 
be also compatible with our finding of /i = 8. The results of Burko and Khanna for 
s — +2 [51] can not be compared with ours due to the different coordinates used there. 

For the first time, to our knowledge, we have verified the LPI splitting (in time) in 
upper modes of s 7^ perturbations, as well as the LPI splitting (in space) of s > 0, 
m = perturbations. In the latter case, the LPIs reach their asymptotic value first at 
a radius R close to the horizon. The LPIs measured by i?+ < R < R show a quahtative 
difference for s = +1 and s — +2 perturbations: in the former case a splitting in time 
is observed very close to it!+. 

In a second series of experiments we investigated the late time decay of non- 
axisymmetric (m 7^ 0) perturbations, mostly of spins s = 0, —2. Focusing on non- 
stationary and compact support initial data (IDl), we considered the late time behavior 
of the field for different values of the black hole spin. For a = 0.9 we find power law 
decays which agree with the literature results for s = and with [44, 46] for s = —2. 

If the black hole rotation approaches the extremal hmit, the late time decays follow 
a different law. Perturbations of a rapidly rotating non-extremal background are at very 
late times still dominated by a weakly damped QNM pattern. In the extremal case, the 
field decay measured from the simulations is consistent with a power law. These 
results agree with the findings of [63, 29], the case s = — 2 is here confirmed for the first 
time by a numerical time-domain experiment. 

In conclusion, the method introduced here is found to be accurate and robust for 
the solution of the TE in the time domain. To the best of our knowledge these results are 
the first numerical investigations of the late time decay at future null infinity for s 7^ 
perturbations on Kerr. We also view these results as preparatory to the development 
of an accurate black-hole-binary test-mass laboratory [72, 73, 74, 62, 75]. Some of the 
results presented may be of interest for the self-force hterature. 

Acknowledgments 

We are indebted to A. Zenginoglu for discussions and pointing out the rescaling in 
Eq. (12) for BL coordinates in an early stage of this work. We thank L. Barack 
and S. Hod for helpful comments on the manuscript. Our work also benefitted from 
discussions with D. Hilditch and A. Weyhausen. This work was supported in part by 
DFG grant SFB/Transregio 7 "Gravitational Wave Astronomy". E.H. wishes to thank 
the council of Off and W. Buchholz. S.B. thanks J.D. Bernuzzi for a careful reading of 
the manuscript. 



Numerical solution of the 2+1 Teukolsky equation, application to late-time decays 33 
Appendix A. Equation coefficients and RT metric 



The coefficients of Eq. (14) are 

_ {R^ -if -M (i?2 _ 1) R^s + 1) + Rh (s cot2(^) - 1) 

R^ 

i?(i? + l) {R^ + lf 
{-a'(i? - l)(i? + If (2M (3i?^ + 5i?^ - 7i?2 - I) - R^ - 6i?' - l) 

+2ia(i? + 1) {R^ + lf s cos{9) 

-2 [4M\R - 1)R{R + If {R\s + 2) + 2R\s + 3) + s) 

+M(i? + 1) {R%7s + 3) - 2i?'\s + 13i?^(s + 1) - AR^s + 2) + i?2(5s - 7) - 2i?s - s - l) 
+2{R - l)i? (i?^s + i?^ + 2R\s + l)+R + s)]} 

C ~ 

{a^ (2i?^ + 5i?2 + 1) 

+2i? [M {R^ - 1) + 3) + 2R^(s + A) + s + l) +2R {R\s + 1) + R^(2s + 3) + s)] } 

cot(^) 



C, 



R 

aR^ + a- 2iRs cot(^) csc(^) 



C'tt — 



2i?(i?+l)(i?2 + l)2 
{a^(-(i? + 1)) (32M^i?^ - (64M2 + 32M + l) i?^ + (32M2 + 32M + 6) i?^ - l) 

W{R+l) {R^ + l)%os(2^) 

-SR [16M^{R - 1)R^{R + if + 8M^R {R^ - R'^ - 3R - l) 

+M [R^ - IIR^ -5R-1) - R{R + 1)] } 

(i?2 _ if (a^ (i?2 _ if + 4^ ^^2 _ ;l) + i?)) 

^^^^ 4i?(i?2 + l)2 

csc^(g) 
~ R 

ja" (i?" - 1)' (2M - 1) - 1) + 2 (4M2i? (i?^ - l)' + M {3R^ + AR + l){R- if - 2R^^ } 
_ 4a(l-2M(i?2-l)) 

__aCR^-lf 
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The coefficients of Eq. (15) are given by 

Cq^ CQ-\-i m §?{C^} - m Q^{C^} - C^^ 
Ce — Ce 
Cee — Cee 

Cr = CR + im Cr^ 

Crr = Crr 

Ct = Ct + i m Ct^ 
Ctr = Ctr 
Ctt ~ Ctt 

Tlie Kerr metric in RT coordinates can be written 

ds" = {-a^+/3RP'^+P^/3'^)dT^+2{PR dT dR+P^ dT dip)+{-fRRdR^+-f0ede^+-f^^d^^+2-fR^ dR dip) , 

(A.l) 

where the 3 + 1 metric functions arc given by 
n {SM\R - l)R + 2M{R + 1) (p2(7? _ 1)2 _ i) _ p^^R _ i)) 

p2(i?- l)3(i?+ 1)2 

_ _AaMRsin^{e) 

p2 (1 - R2) 

4 (4Mi?2 - (4M + 1)R - 1) {16iVPR^{R - 1) + AMR{R + 1) {p'^{R - if - 1) + p^R - 1)^) 



IfRR 



p\R-i)^R + iy 

2 



lee = p 

_ S^sin^(g) 



_ 2a sin^(g) {16M^{R - 1)R^ - AM{R + 1)R + p^ {R^-R^ + R- 1)) 

'^'^^ ~ p^iR-ifiR + iy 



a 



p\R - If {R' + ly (a^sin^(^) - - 

sin^(^) (l6M^(i? - l)i?^ - AM(R + 1)R + p'^ {R^ - R^ + R - l))^ 
+{R -1){R+ 1)S^ (4Mi?2 - (4M + l)i? - l) 

{16M^R\R-1)+AMR{R+1) {p^{R-lf - l)+p^{R-lf)]~^y 



and the shorthands 



1 -i?2 

A{R) = r{Rf -2Mr{R) + 

E{R, e) = ^ (r(i?)2 + a2)2 - a2A(i?) sin(e)^ 
p{R, 9) = Vr(i?)2 + a2cos(^)2 
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are used. 

Appendix B. Green's function calculation 

In this appendix we review the analytical calculation of [45, 44, 46], and specify it for 
the case of IDl pointing out explicitly a special case relative to m = 0, s 7^ and pure 
multipole initial data I' — Iq + 1 that was overlooked in the final parametric formulas 
of [46]. To avoid confusion we will follow the naming conventions of [46], but keep I' 
as index of the pure multipole initial data. The calculation is in BL coordinates, and 
employs the radial tortoise coordinate y defined by dy = (r^ + o?')/Adr. The main 
variable \1/ refers to the 2+1 decomposed fields rescaled by A~*/^(r^ + a^)~^/^ (compare 
with Eq. (12).) 

The solution to the 2+1 TE wave equation for t > can be expressed in terms of 
the retarded Green's function G{y, 9, y', 9', t) and the initial data as (see Eq. (6) in [46]) 

*(y, 9, t)^2nj {B^{y\ 9') [G(y, 9, y\ 9\ t)^^{y\ 9\ 0) + G,(y, 9, y' , 9', t)^(y', 9', 0)] 

(B.l) 

+B2{y\9')G{y,9,y\9',t)^{y',9\Q)} sm9'd9'dy' , 

where Bi{y,9) and B2{y,9) are background dependent complex functions (Eq. (4-5) 
in [46]), in particular Bi{y,9) — 1 — bi{y) sin^^. 

For non-stationary, pure multipole IDl '^{y,9,0) = with ^^(^,^,0) oc ^Yirm(9) 
only the first of the three terms in Eq. B.l has to be evaluated, 

^{y,9,t) = 27T J ^ B,{y',9')G{y,9,y',9',t)^t{y',9',0) sin 9' d9'dy' . (B.2) 

The Green's function G{y, 9, y' , 9' , t) can be expressed in terms of its Fourier transform 
and expanded in the spin weighted spheroidal harmonics ^Sim{9, auu) basis, 

G{y,9,y',9',t) = {27t)-^ j Y.^i{v.0,y\9' ,ujySU9,aujrSU9\auj)e-'^'d^ . 

-oo+ic ^='0 

(B.3) 

The Green's functions Gi{y,9,y' ,9' admit an analytic solution in terms of two 
linearly independent solutions of the radial Teukolsky equation, sec Eq. (18) of [45]. 
Given Gi{y,9,y' ,9' ,uj), to obtain the Green's function it remains to plug it back 
into Eq. B.3 and evaluate the integral. As discussed in [76] the integration contour 
is chosen in the lower half of the complex u plane. The integral consists of three 
distinct contributions: (i) an integral along the semicircle contributing to the early 
time solution, (n) the residues (from the poles of Gi{y,9,y',9',u)) in the lower half of 
the plane, corresponding to the QNM contribution; (iii) an integral along the branch 
cut Gi{y,9,y' ,9' ,uj) on the negative imaginary uj axis, contributing to the late time 
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solution. In the following wc consider only the latter and write the Green's function as 
G{y,e,y',e',t)^G^{y,e,y',e',t). 

The calculation of the tail contribution to the Green's function, and thus of the 
solution in Eq. B.l, proceeds by considering a small u approximation, which corresponds 
to selecting the large r contribution of the (effective) scattering potential. Gonsidering 
Gi{y,9,y' ,6' ,u!) in the small-u; hmit the expression for G"^ for observers at time like 
infinity reads (Eq.(12) of [46]) 

G^{y,e,y',e',t) = J2Kis{yyy^' J 'Si,n{e,aoorSi^{e',auj)uj^'+'e-'^' dw (B.4) 

l=lo 



where Kis are constants solely depending on / and s. 

Next, appropriate approximations for ^Skmi^, a,cj) have to be introduced. In general, 
the spheroidal spherical harmonics are defined by the eigenvalue problem 

where L is a differential operator which splits into a a;-independent part with 9- 
derivatives (Lq) and Li — {auY cos^ 9 — 2{au!)scos9, the a;-dependent part. If (aw) 
is small, one can view Li as a perturbation of the operator Lq and compute a 
perturbative expansion in {acu) [77]. The unperturbed eigenfunctions are the spin 
weighted spherical harmonics, which are the eigenfunctions in the non-rotating case, 
i.e. '^Sim = ^Yim + C(aa;). The expansion then reads 

where {sim\F(9)\slm) = J dfl ^Y*^^YimF{9), and AA^^-* is the difference between the 
Z-th and i-th unperturbed eigenvalues (independent of 9 and uj). The computation 
of the expansion coefficients requires the evaluation of the integrals {sim\ cos 9\slm) 
and (simj cos^ ^jsZm), the former are non-vanishing for i = 1,1 ± 1, the latter for 
i — 1,1 ±1,1 ±2. The expansion in Eq. (B.6) can be written as 

'Sirn = Y,Cu iauf-'^ 'Y,^ , (B.7) 

i=lo 

where the coefficients Gu depend on {aou), but are at leading order independent of (auj) 
for (au) <^ 1. 

Substituting the unperturbed eigenfunctions in Eq. (B.4), i.e. ^Sim — > *5^m, leads 

to 

—i oo 

oo „ 

G^(y,9,y',9',t) = J2Kis{yy'y+''YU0rYrm{0') J e— * da; oc r^^'^^) , 

(B.8) 
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which reproduces the Price law with /j, — 21 + 3. In order to describe mode couphng 
and obtain the actual tails on Kerr, one needs to consider the expansion in Eq. (B.6). 
Renaming some indexes, dropping the s and m indexes, and rearranging some terms, 
Eq. (B.4) reads 

—i oo 

oo oo oo „ 

= E E E J dcu e-'-' ^2^+2 C,i {aujf-'\ Cu {owf-'^ , 

1=Iq i=lo k=lo Q 

(B.9) 

from which one can read off the Z-mode contribution Gf defined by = Gf. 

We define now an effective Green's function for the Z-mode that (i) takes 
into account the non-vanishing terms in the sum over i to the (outer) integral 
/ de' sme'Bi{e'yYi,m{6'), and (ii) restricts to the lowest powers of (auj) that contribute 
to the sum over k. Step (i) requires the evaluation of the integrals like {sim \ sin^ 6'\sl'm). 
In particular, for m 7^ one has that the non-vanishing terms of {sim\ sin^ ^^'|s/'m) are 
those given by i = /' ± 1, /' ± 2. For m = one has that the non- vanishing terms of 
(siO| sin^ 9'\sl'Q) are those given by i = ± 2. The calculation at step (i) gives 

^ —ioo 
00 „ 

Gf'^^Y.^'^^^yy')''^' j dwe-'^' u^^^^Cu{owt-'\Yi{e) (B.IO) 

k=lo Q 

{[l-h,{y')]Cu'{auot-'%r{e') 

- h{y') [ c,i,_2iacof-^'-'^\Y;_,{e') + c,i,_,iawf-^'-'^\Y;_,ie') 

Note that of the terms in the curly brackets the Ckv±i do not exist for m = and 
the Ckv-2 and Ckv-i do not exist, respectively, for /' = /o,/o + 1 and /' = /q- Step (ii) 
requires power counting for given combinations of (/, /'). The freedom of the powers in 
the curly brackets requires to analyze distinctly the three cases I' — Iq, I' — Iq + 1 and 
^' > ^0 + 2. In each of these checking the first terms of the k-sum reveals that the lowest 
power is always obtained by the k — Iq term (however k > Iq terms may amount to the 
same power). Thus the lowest power of (ao;) is 

'2/o + 2 + |/o-/'| + |/o-/| , ^' = /o , 

=<( 2/0 + 2 + |/o- (/'- 1)1 + |/o-/| ,/' = /o + l, (B.ll) 

2k + 2 + |/o - (/' - 2)1 + |/o - /| , /' > /o + 1 . 
The lowest power of (ao;) for m = is instead given by 

2Zo + 2 + |Zo-/'| + |/o-^| , ^' = ^0 , 

ni^ {2h + 2 + \h-l'\ + \h-l\ ,r = Zo + l, (B.12) 

^2/0 + 2 + |/o - (/' - 2)1 + |/o - /| , r > /o + 1 , 

as a consequence of the vanishing integrals (siO| sin^ ^^'|s/'0) for i = l'±l. It is important 
to note that for I' > Iq + 2, other terms in the sum than the k = Iq can occur with the 
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power Til. This is because, while the u '^'^ contribution will always be greater than 
1^2/0+2 ^Qj, k > Iq, the (jj powers that occur with the C coefficients can be smaller for 
k > lo than for k = Iq. 

Finally, the power of the tails for each /-mode can be calculated from the integrals 
J da; e~*'^*C(;"' which give with /x; = ui + 1. The different cases above are actually 
summarized in the parametric formula in Eq. (15) of [46]. However, the parametrization 
there seems to miss the special case I' — lo + l and m — 0. 

The calculation at null infinity is similar to the one above. In this case different 
approximations for the Fourier transform Gi{y, 9, y' , 6' , u) of Eq. (B.3) have to be made 
but the remaining procedure is unchanged. In particular using again Eq. (B.7) for the 
approximation of the ^Sim and that {sim\ sin^ 6'\sl'm) = unless i — I', I' ±1,1' ±2 one 
finds the effective Green's function (cf. [46]) of the /-mode 

Gf'^ = J2 Kks y~' j dw e-'^^'-y^ uj''-'+^ Cki (acu)!^"'! YiiO) (B.13) 

k=lQ Q 

+ Ck,+,{auf-^'+'^\Y;^,{e') + c,,+2{aujf-^'^'^\Y;^,{e') ] } , 

where Kks is a constant depending on k and s. This Green's function provides again all 
information to find the decay rates for any pair of (/, /') by going through the sum over 
k and picking the lowest power of u. While at the lowest power is always obtained 
by the k = Iq term one finds that this is not true at J^+. Still in [46] Hod found an 
appropriate way to parametrize the resulting decay rates but again the parametrization 
contains an error for /' = /q + 1 with m = 0. In this case the Cfc/'_i(aa;)'^~*^''~^'"y^?_]^(6'') 
term in the curly brackets vanishes so that for / = Iq Hod's prediction is an integer 1 
below the actual value. Higher /-modes of the I' — Iq + 1, m — case are predicted 
correctly because in these cases the k — Iq + 1 term gives the same power of ou without 
usmg the Cw/_i(aa;)l'=-(''-^)lyjr_i(e') term of the curly brackets. 
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